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This study was a sensitivity analysis of a computer 
code (based upon CHAMP ION/ 2/E/FIX ) , used at the Naval Post- 
graduate School mainly to model the flow within Solid Fuel 
Ramjets. The purpose of the study was to make the code 
"simpler" and "more accurate" by reducing the required 
CPUs time and/or improving the accuracy of predictions . 
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INTRODUCTION 



I . 

During the last fifteen years, there have been many 
advancements in the development of highly effective tacti- 
cal weapons for use at intermediate range and with high 
average flight speed. One of the simplest engines/motors 
for such use is the Solid Fuel Ram Jet (SFRJ) . 

The SFRJ most often consists of a solid fuel grain which 
provides the walls for the combustion grain [Ref. 1] . A 
sudden expansion at the air inlet of the combustor can be 
used to provide flame stabilization by inducing recircu- 
lation into the air flow (Fig. 1.1). . Combustion can be sus- 
tained throughout the grain provided that the step height 
(h) is sufficiently large. However, the larger the step 
height, the larger the losses of the inlet stagnation pressure. 

Due to the sudden expansion, at least two distinct flow 
fields are generated within the fuel grain (Fig. 1.2). In 
the first field, the recirculation zone, the flow is highly 
turbulent and usually fuel rich. This hot gas region pro- 
vides the energy necessary to sustain the combustion proc- 
ess which must occur further downstream. Downstream of the 
flow reattachment, a turbulent boundary layer develops and 
the combustion is diffusion controlled. A diffusion flame 
emanates from the recirculation zone and remains within 
the developing turbulent boundary layer between the 
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Figure 1.1. Schematic Diagram of Solid Fuel Ramjet 




Figure 1.2. Effects of Sudden Enlargements on a SFRJ 

fuel-rich zone near the wall and the oxygen-rich zone near 
the central core. 

There have been continuing research efforts concerning 
the behavior of the SFRJ at the Naval Postgraduate School. 
Both mathematical and experimental efforts have taken place 
to determine the effeclis of design as well as operational 
variables of the system to be obtained. 
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II. DESCRIPTION OF THE COMPUTER MODEL 



A. INTRODUCTION 

At the NPS there have been two basic computer codes 
used to model the flow within solid fuel ramjets. The 
first one was based on stream function ( 4 ;) and vorticity (o)) 
[Ref. 3]. It predicted with reasonable accuracy the velocity 
and temperature distributions but did not predict pressure 
distributions accurately and was not easily extended to 
complex geometries or to three-dimensional flows. The second 
computer model is based on the "primitive” variables, pres- 
sure (p) and velocities (u,v) [Ref. 2J . It has also been 
used to predict the effects of fuel properties on the SFRJ 
performance and to evaluate the effects of different 
geometries as well as operating conditions. 

B. OVERVIEW OF THE PRIMITIVE VARIABLE PROGRAM 

1 . Assumptions 

The CHAMP ION/ 2/E/FIX computer program developed by 
Pun and Spalding [Ref. 2J has been used as the basis for the 
primitive variable model. The flow has been assumed to be 
steady (solid fuel grain dissipation is negligible! , sub- 
sonic, recirculating, axisymmetric and to have constant 
specific heats. In order to calculate the effective viscosity, 
a modified Jones-Launder [Ref. 4J two-parameter turbulence 
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model has been used. It uses five empirical constants 
(Table II) and requires that two additional variables 
(turbulent kinetic energy (k) and turbulence dissipation 
rate (e)) be evaluated. The effective viscosity was calcu 
lated using the formulas 



''eff = '‘lam * "t 



( 2 . 1 ) 



where 



Wt = Cp p K /e 



( 2 . 2 ) 



TABLE I 



k-£ Turbulence Model Empirical Constants 



k,eff 



e , ef f 



1.43 



1.92 



.09 



1.0 



1.3 



The combustion was considered to be mixing- limited 
(infinitely fast kinetics) with a simple one-step chemical 
reaction of the form: 

1 Kg fuel + i Kg, oxidizer (l+i)kg products (2.3). 

Four species were concerned: oxygen, nitrogen, fuel 

and products. Because of the infinitely fast kinetics 



17 



TABLE II 



Governing Equation Parameters 




Cl = 1.43, O 2 = 1.92, \ = 1' 
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assumed, no oxygen could exist at the fuel surface and at 
the surface was assumed to be isothermal . The turbulent 
Prandtl and Schmidt numbers were taken to be equal to unity 
and therefore the turbulent Lewis number and laminar Prandtl 
number were also taken to be equal to unity. 

C. GOVERNING EQUATIONS 

With the above assumptions, the governing equations for 
axi symmetrical flows can be written as follows: 



|j(pu<|,) + i |j(prv4.) - 



— ( r 

8x^' (}) 



- k 1 



8x 



(rr 

8r ^ (p 



d <p 

Jr 



) + S . = 0 (2.4) 

9 



where stands for the independent variables u, v, k 



h, ^ first two terms are the 

tion terms", the next two the "diffusion terms". F 



<P 



f 

"convec- 
is the 



effective exchange coefficient for turbulent flow (Table I) 
and S^ is the "source term" (Table II) . The stagnation 
enthalpy for unity Prandtl and Schmidt numbers is 



h = h + C'u^+v^)/2 + K 



C2.5). 



where for nonreacting flows 



h 



C T 
P 



C2.6) 



and for reacting flows 
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h 



(2.7) 



= m ^+c(T-T^) 
ox 1 p ref 



Temperatiore has been calculated using Eqs. (2.5) 
and (2.7). Density has been calculated using the perfect 
gas law 



= 

RT 



( 2 . 8 ) 



Since the quantities h (stagnation enthalpy) , 



and Y = - m /i 

fu ox^ 



have identical governing differential equations and, in 
appropriate dimensionless form, identical boundary condi- 
tions, a representative equation had to be solved only once. 
The dimensionless enthalpy was selected to be 



H = 




C2.9) 



where h^^ is the stagnation enthalpy deep in the fuel grain 
and h^^ is the stagnation enthalpy of the air at the inlet 
dump plane. The corresponding dimensionless forms for the 
other quantities were: 



and 




X- 

2 in 









- m, 



2 in 



N 



2f g 



C2.10). 
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X 




( 2 . 11 ) 



m 



fu 



X 



m 



0 



for X > 0 



( 2 . 12 ) 



m 



fu 



0 



m 




for X < 0 



ox 



1 . Boundary Conditions and Solution Procedure 

Fixed boundary conditions (inlet velocity profile, 
et al.) were specified at the desired or experimentally 
determined values. Specified gradient boundary conditions 
were handled by setting the appropriate convection/diffusion 
coefficients to zero in the finite difference equation 
("breaking the link") and then entering the appropriate 
gradient through linearized "false" source terms. The. 
geometry as well as the appropriate boundary conditions are 
summarized in Figure (2.1) . 



this was not a computer program limitation. The turbulence 
kinetic energy was assumed to be uniform, with a value which 
corresponded to the approximate turbulence intensity of the 
inlet flow. Radial and axial gradients were set equal to 
zero on the centerline and at the exit, respectively. All 
the nonreacting solid boundaries were considered adiabatic 
with both u and v velocity components equal to zero (no-slip 
condition) . 



•I 



Plug flow" was assumed at the air inlet, although 



21 



o 



o ^ 

*» 




o 

« 



ek 



o 

n 

? 

II 



<• 

o 







f o- 
v*> 



!• a 



O 

O ” 

U ^*e 

? (?‘C. 



V 















o 




o 








M 

0^ 

r 


'S 

o 


O 

? 
5C * 


o 




o 


o ' 


M 


1* 


• 

n 


• l» 





w 

0) 

jC 

■u 

M 

O 

CO 

C 

o 

•H 

■P 

•H 

'O 

c 

o 

u 

>1 

u 

(d 

'O 

c 

3 

o 

(X) 

ip 

0 

>1 

u 

+J 

0) 

e 

0 

0) 

u 



04 

0) 

3 

•H 

b 



22 



For simplicity, the boundary layer was assumed to be 
composed of a laminar sublayer which is then developed to 
turbulent without assuming a "buffer zone". For the 
laminar sublayer. 



= y. 



(2.13) 



where 



u 



u 



/t /p 

w 



and 



+ 



’'Vpy, 



V 



and for the rest of the layer 



u 



iln(Eyp) 

K 



(2.14) 



For E = 9.0 and k = 0.4 (values found from experimental 
data) 



u"^ = 5.5 + 2.5 iln y'^ (2.151 

P P 

Solving the above two equations, it is found that 

transition from laminar to turbulent flow occurs at y"*" = 11.5. 

The value of y'*’ at each near-wall node is found from the 

P 

relation: 





^lam 



w 

P 



1/2 



(.2.16) 
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For y < 



11.5 



T 

W 



y, * u 
lam p 

6 



(2.17) 



and for 




> 



11.5 



T 

W 




(2.18) 



The value of x was assumed to be uniform from the wall to 

w 



the near-wall grid point. So, 



+ 



y 



p 




p 



1/2 6 
K ' 

P ‘'lam 



(2.19) 



Then, the shear stress is calculated using the formula 
[Ref. 3] 



w 



1/2 

D p 



u 



= p c 



1/4 ^i/2,::£, 



-^) 
p + 

p 



E D c^/^ 

Kpu in{ 2 ) 

P P l^lam 



( 2 . 20 ) 



where, as noted previously. 



u 



In (Eyp) 



Due to the steep gradients of the properties in 
turbulent flows near the solid boundaries, the source terms 
for K and e at the near-wall nodes were expressed in terms 
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of the wall shear stress. Also, provides the boundary 
condition for the u and v equations. 

The turbulence dissipation rate at a near-wall 

node p is given by the relation [Ref, 3] 

Sp = = <p'^V2.4 36 (2.21) 

where the length scale is presumed to be proportional to the 
incremental distance from the wall, 

£ = c^^'^ K 6 (2.22) 

It was found that when using a sudden expansion 

geometry in reacting flows the near-wall dissipation rate 

3/2 

had to be increased on the step face (e = k ' /0.4(S) and the 

P P 

grid spacing adjacent to the grain had to be fine (,y^ < 11,5) 
in order to obtain a temperature distribution in qualitative 
and quantitative agreement with experimental data [Ref, IJ . 
Equation C2.17) implies that the wall shear stress is calcu- 
lated assuming a linear velocity profile when y^ < 11.5. 
Therefore, a near-wall grid point can lie within the laminar 
siablayer but the source terms for k and e imply that 
is much greater than unity. This precludes y^ from being 
significantly less than 11,5. 

For reacting flows, the boundary conditions for the 
dimensionless properties (Eqs. (2.9) through (2.12)) were zero 
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at the inlet and unity "deep" in the fuel grain. These 
properties were considered to have zero gradients on non- 
reacting surfaces. 



The assumptions made for reacting flows (unity 



Prandtl and Schmidt numbers, simple chemical reaction, con- 
stant specific heat and stagnation enthalpy defined in Eq. 
(2.5)) result in a general boundary condition for all 
"conserved" properties (<{) ) on a surface which has mass 
transfer given by the relation [Ref. 2] 



where (|)^ may represent h,m^ or x = m^^ - m^^/i. 

2 

A mass transfer conductance (g) is often defined 
from the. relation 



where 4)^ is defined as the free-stream value of the property. 

Substituting Eq. (2.24) in Eq. (2.23), the mass 
transfer from the blowing wall is 







(2.23) 



^ ^ <j) dr bw 





(2.24) 




^^^cp ^cbw^ ^^c,bw ^c,fg 



C2.25) 



where 
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BP 



(2.26) 



^cp ^C/bw 
^c,bw ^c,fg 

represents the mass transfer parameter (or blowing parameter) . 
The value of g can be approximated using 



g = (pu) St 



(2.27) 



From Reynolds' analogy with unity Prandtl number, the 
Stanton number is 



ST = 




w 



/ 2 , 

(pu ) 



(2.28) 



where is the local friction coefficient. 



Thus , 



g 




(2.29) 



Using the Couette flow approximation for the boundary- layer 
behavior with mass transfer, 



g = g* £n(l + BP) /BP 



C2.30) 



where 



* 

g 



lim (g) 
BP-^0 



(2.31) 
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In this case, BP has been evaluated from the solution 



of the energy equation using 

h - h, 

BP = (2.32) 

bw fg 

The wall shear stress has been evaluated using 
Eq. (2.17) or Eq. (2.20) and modified with Eq. (2.30) 



T, = T £n(l + BP) /BP 
bw w 



(.2.33) 



where x is the wall shear stress without wall mass addition, 
w 

The mass transfer conductance (g) was evaluated 
using Eq. (2.29). The wall heat flux (q") on all solid 
isothermal boundaries was evaluated using Reynolds ' analogy 






(2.34) 



Since the blowing rates were small for the solid fuel 
ramjet (typically BP < 2.0), k and e were evaluated using 
Eq. (.2.4) and the terms presented in Table II, which incor- 
porate the empirical constants of Table I. 

Five variables, namely u, v, k, e and H Cor) h were 
solved using Eq. (.2.4) in finite difference form. The line- 
by-line iterative procedure employed upwind differencing 
and underrelaxation in order to promote convergence. On 
each radial line, the mass flow rate was calculated using the 
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local density. The "error" in mass flow (compared with the 
summation of "mass- in" at the upstream boundaries) was used 
to uniformly adjust the axial velocity over the entire line. 
Then, the pressure at all downstream locations was adjusted 
to approximately correct for the momentum imbalance created 
by the uniform axial velocity and a "pressure correction" 
equation was solved for each cell on the line. Finally, 
local cell velocity (axial and radial) and pressure were 
adjusted to satisfy cell-wise continuity. 
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III. ANALYSIS 



A. INTRODUCTION 

As previously discussed, the purpose of this study was 
two-fold; first, to reduce the CPU time required to run the 
program, without considerably affecting the accuracy of 
results and second, to make the program output better fit the 
experimental measurements without considerably increasing the 
CPU time. The methods used in this optimization procedure 
were the following; 

a. Monitor the effects of the number of traverses in the 
tridiagonal matrix (TDMA) . 

b. Monitor the effects of the nxmber of iterations. 

c. Monitor the effects of the grid spacing close to the 
solid fuel grain wall, without changing the number of 
grid lines, in either the x or the r directions 
(constant number of nodes) . 

d. Monitor the effects of changes in the number of grid 
lines in both the x and r directions Cthe number of 
nodes was changed) . 

B. EFFECT OF TDMA TRAVERSES 

The number of TDMA traverses in the main program is denoted 
by the parameter NTDMA. If there is a predominant direction 
of flow, this is arranged to be in the x-direction. Then, 

TDMA traverses are made on constant x lines, in order to make 
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the adjustments necessary. Generally, the greater the number 
of TDMA traverses, the better the accuracy obtained but the 
more CPU time is consumed. 



C. EFFECT OF NUMBER OF ITERATIONS 

The number of iterations in the main program is denoted 
by the parameter LSWEEP. A high number of iterations (or 
sweeps) is needed in order to have small residual values in 
the variables used. Here again the same difficulties as with 
the nTomber of TDMA traverses appeared. The larger the value 
of LSWEEP, the better the accuracy but the greater the CPU 
time required. It is expected that for a high number of 
iterations, the increase in accuracy will be small whereas 
the CPU time will be excessive. 

D. EFFECT OF GRID SPACING CLOSE TO THE SOLID FUEL GRAIN WALL 
In the original program, the grid lines spacing close to 

the physical boundaries was small. For small grid spacing 

close to the solid fuel grain wall, the non-dimensional length 

y^ at near wall grid points was much less than 11.5 and the 
P 

flow was assxmied to be laminar. The velocity profile in this 
case is linear, given by 





^lam 



(.3.1) 



For dense packing of grid lines, the flow may be con- 
sidered as laminar for more than the first node close to the 
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wall. In this case, the velocity profile will consist of 
several straight lines before turbulence will be encountered, 
as shown on Fig. (3.1). This could result in inaccurate 
predictions. The spacing was changed such that the first 
constant-r grid line was spaced far enough from the wall 
in order to have > 11.5. In this case, the velocity pro- 
file will be a logarithmic curve. If this increased spacing 
resulted in good predictions, then the number of constant-r 
grid lines could be reduced, thus reducing the total CPU 
time. 

E. EFFECT OF CHANGES IN THE NUMBER OF GRID LINES 

Generally, the greater the number of grid lines in both r 
and x-directions, the better the accuracy expected but also 
the greater the storage as well as the CPU time required to 
execute the program. Since it is desirable (but not neces- 
sary) to have a maximum ratio of 1.5 for distances between 
adjacent grid lines in the same direction and a maximum length 
to width ratio of 10.0 for any cell, it is preferable that 
both the number of constant-x and constant-r grid lines be 
changed (denoted by NX and NY respectively in the main program) . 
In such a case, if the grid spacing close to the solid fuel 
grain wall could be increased — as discussed above — the number 
of lines in both directions could be reduced without affecting 
the accuracy of predictions. This would result in reduced 
requirements for storage and CPU time. 
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F igure 3.1. 



velocity Profile Close to Solid Fuel Grain 
Wall 



33 



IV. RESULTS 



A. INTRODUCTION 

All the options discussed in the previous chapter have 
been tried in the modified CHAMP ION/ 2 /E/FIX computer program. 
The original program had been modified in the past at the 
Naval Postgraduate School to account for an aft-mixing 
chamber [Ref. 2 ], and for compressibility and radiation 
effects [Ref. 6] . This program calculates the flow in a 
SFRJ, whose dimensions are shown on Fig. (2.1). The input 
values used in the program are shown in Appendix B, which 
is a copy of the BLOCK DATA subroutine. 

The results of interest are plotted and discussed in 
the following sections of this chapter. The plots mainly 
refer to two x-constant grid lines after the first and second 
sudden enlargements at x = 202 mm and x = 312 mm, respectively. 

B. EFFECT OF NUMBER OF TRAVERSES 

The original program has been modified to run with three 
different values of TDMA traverses: NTDMA = 3, NTDMA = 5 

(originally used value) and NTDMA = 7. Nothing else had 
been changed. An attempt to run the program with NTDMA 
greater than 7 resulted in a repeated series of overflows 
after a relatively high number of iterations had taken place 
(LSWEEP > 20) and the program execution was terminated. 



34 



1 . Temperature Distribution 

The temperature profile variations were negligible 
in this case, as shown on Figs. (4.1) and (4.2). Thus, by using 
a smaller number of traverses, the accuracy of temperature 
distribution prediction was not significantly affected. 

2 . Effective Viscosity Distribution 

Here again, the effect of variations of effective 
viscosity, as shown on Figs. (4.3) and (4.4), were negligible. 

3 . Axial Velocity Distribution 

The axial velocity distribution after the first and 
second sudden enlargements are shown in Figs. C4.5) and 
(4.6), respectively. Once more, the number of TDMA traverses 
did not affect the accuracy. 



C. EFFECT OF NUMBER OF ITERATIONS 

The program was run for three different numbers of 
iterations; for LSWEEP = 100, LSWEEP = 200 (original value! 
and LSWEEP = 400, without any other changes. 

1. y^ Distribution 

The values of y^ after the first sudden enlargement 
were essentially unchanged. Some differences appeared after 
the second sudden enlargement but these results did not seem 
to depend upon the number of iterations; the peak value of 
y^ for LSWEEP = 100 was between the peak values for LSWEEP = 200 
and LSWEEP =400. In any case, the value of y^ defies 
accurate prediction because of recirculating flow; the relatively 
high value of y^ was an indication of high turbulence. Thus, 
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changing the Number of TDMA 
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Enlargement, Changing the Number of TDMA 
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Figure 4.6. Axial Velocity Distribution After the Second 

Sudden Enlargement, Changing the Number of 
TDMA Traverses 
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a change in the number of iterations did not significantly 

affect the y**" distribution, 

P 

2 , Temperature Distribution 

The differences in the temperature distributions after 
the first sudden enlargement were negligible, as shown in 
Figure (4,8). The peak temperature, 2700°K, which was an 
indication of flame location, moved almost to the mid-radius 
point. Then, a relatively steep reduction of temperatxire 
occurred, and at the centerline its value was less than 
2000®K. 

After the second sudden enlargement, the temperature 
distribution differences were also negligible, as shown in 
Figure (4,9). The peak temperature now occurred at the center- 
line of the SFRJ, which meant that the enlargement effectively 
promoted the mixing of oxygen and unburned fuel. 

3. Effective Viscosity Distribution 

Figure (4.10) shows that the variations in effective 
viscosity distributions after the first sudden enlargement 
were almost negligible. All three curves coincided close 
to the solid fuel grain wall and then had a small separation 
as they approached the centerline. The same effect can be seen 
in Figure (4.11). These small differences in close to 

the centerline did not really affect the calculations since, 
as will be seen later, the axial velocity distribution close 
to the centerline was amost uniform. As a result, the shear 
stress value was also low and was not affected by the variations 
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Figure 4.8. Temperature Distribution After the First Sudden 

Enlargement, Changing the Number of Sweeps 




(W) 3didy 



45 



Figure 4.9. Temperature Distribution After the Second Sudden 

Enlargement, Changing the Number of Sweeps 
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Figure 4,11. Effective Viscosity Distribution After the 

Su(3(3en Enlargement, Changing the Number of 



From the above, it is seen that a reduction in 



of Pgff. 

the number of sweeps to one half of the originally specified 
value did not significantly affect the effective viscosity 
distribution. 

4 . Axial Velocity Distribution 

As for temperature and effective viscosity distribu- 
tions, the variations in the distribution of axial velocity 
were also negligible for this range of iterations. Close to 
the solid fuel grain wall, the values of axial velocity were 
the same in all three cases. Close the centerline there was 
a separation, with higher values attained for higher numbers 
of iterations for both of the areas after the first and the 
second enlargement Cshown on Figures (.4.12) and (4.13), 
respectively) . 

Here, it was evident that increasing the number of 
iterations increased the accuracy attained; although, 
for this range of sweeps, the variation in the predictions 
did not exceed 3%. 

D. EFFECT OF GRID SPACING CLOSE TO THE WALL 
1 . y^ Distribution 

The grid spacing was increased close to the solid 
fuel grain wall, such that turbulent conditions would exist 
at the near wall grid points all along the combustion length. 
The values of y^ have been plotted in Figure (4.14),. Values 
higher than 11.5 indicate turbulent flow, whereas lower values 
indicate laminar flow. 
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Figure 4.12. Axial Velocity Distribution After the First Sudden 

Enlargement, Changing the Number of Sweeps 
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Figure 4,13. Axial Velocity Distribution After the Second Sudden 

Enlargement, Changing the Number of Sweeps 
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Fiaure 4.14. Distribution, Changing the Grid Spacing 

cEose to the Solid Fuel Grain Wall 



As shown on Figure (4.14 ) , increasing the grid 

spacing close to the wall resulted in the values of being 

generally increased along the region from the first to the 

second sudden expansion. After the second sudden expansion, 

the values of y^ that corresponded to increased grid spac- 

P 

ing were generally lower, thus indicating less turbulence. 
However, the values were greater than 11.5 and the flow 
was handled as turbulent from the second expansion to the 
exit. 

2 . Temperature Distribution 

The temperature distribution was dramatically 
affected by the grid spacing close to the solid fuel grain 
wall (i.e., the distance from the first node to the wall) . 
Figure (4.15) shows that the original program predicted a 
peak temperature rather far from the wall. On the other 
hand, the program with increased spacing from the wall pre- 
dicted peak temperat-ures close to the wall (the larger the 
grid spacing, the closer the distance of peak temperature to 
the wall). . A peak temperature farther from the wall means 
a thicker boundary layer and, in turn, a more fuel rich 
mixture ratio . 

Figure (4.16) shows that the original program pre- 
dicted complete mixing in the aft mixing chamber . Peak 
temperature had been predicted to be at the centerline of 
the SFRJ, whereas with increased spacing the peak temperature 
occurred just behind the step wall. 
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Figure 4.15. Temperature Distribution After the First Sudden 

Enlargement, Changing the Grid Spacing Close 
to the Wall 
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3 . Effective Viscosity Distribution 



As shown on Figure (4.17), the effective viscosity 
distribution close to the well was not greatly affected by 
the spacing of the grid lines. However, close to the center- 
line the effective viscosity value increased with increasing 
spacing, while the value of 1 -Iq££ close to the wall was slightly 
greater in the original program. 

The same happened after the second expansion, as shown 
on Figure (4.18), although the curves obtained were now 
somewhat different. Again, the effective viscosity from the 
original program was greater close to the wall, while, as 
the centerline was approached, this value dropped again below 
the values predicted for increased grid spacing close to the 
wall . 

4 . Axial Velocity Distribution 

As shown on Figures (4.19) and (4.20), the original 
program predicted much higher velocities than the programs 
run with increased grid spacing close to the wall. It was 
also noted that the velocity profile of the original program 
was typical of laminar flow, where the ones for increased 
spacing were typical of turbulent flows (almost uniform 
velocity close to the core) . 

In this case too, the same general effects in the 
distribution could be noted as in those for effective 
viscosity. The original program predicted higher values than 
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Figure 4.17. Effective Viscosity Distribution After the 

First Sudden Enlargement, Changing the Grid 
Spacing Close to the Wall 
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Figure 4.18. Effective Viscosity Distribution After the 

Second Sudden Enlargement, Changing the 
Grid Spacing Close to the Vv’all 
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Figure 4.19. Axial Velocity Distribution After the First Sudden 

Enlargement, Changing the Grid Spacing Close to the Wal 
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Figure 4.20.. Axial Velocity Distribution After the Second 

Sudden Enlargement, Changing the Grid 
Spacing Close to the Wall 



the programs with increased spacing and the .velocity pro- 
files corresponded to laminar and turbulent flows, respectively. 

These large variations in the u-velocity distributions 
shown on Figures (4,19) and (4.20) were not surprising con- 
sidering the variations in the temperature distributions. 

Since, for increased grid spacing, the temperature was lower 
in the core, the velocity was also lower than that obtained 
in the original program. 

E. EFFECT OF A WIDER GRID SPACING 

Since this case had similarities with the previous one, 

only one case was run; with increased grid spacing in both 

X and r directions. The y^ vs. x distribution has been 

P 

plotted, as shown on Figure (4.21) . It was found that this 

case was not of particular interest. Along the solid fuel 

grain wall y^ was slightly increased. This in turn meant 
P 

that predictions would not be as accurate as those obtained 
increasing the grid spacing close to the wall only. So, no 
further analysis took place for this case. 
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V. CONCLUSIONS 



From the Analysis of Chapter IV it can be stated that; 

a. The computer program was not sensitive to changes 
in the number of TDMA traverses, NTDMA. Thus, a lower 
value, such as NETDMA = 3 can be used without reducing the 
accuracy of predictions . 

b. As expected, the higher the number of sweeps, LSWEEP, 
the greater the accuracy. However, if the number of iterations 
exceeded 200, the increase in accuracy was negligible. For 
any number of sweeps between 100 and 200, the accuracy of 
prediction was found to be very satisfactory. 

c. The program was very sensitive in the grid spacing 

close to the solid fuel grain wall . For increased grid spacing 

» 

in such a way as to have y^ > 11.5 for most of the nodes 
adjacent to the solid fuel grain wall, the results of the 
program fit experimental data of Well [Ref. 5], as shown on 
Figure (5.11 . Also, this increased spacing close to the 
wall permitted the' reduction of the number of grid lines 
in the r and, consequently, the x-directions , with consider- 
able savings in preparation as well as CPU time. 
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with Experimental Data 



VI . RECOMMENDATIONS 



In view of the results of Chapters IV and V, the follow- 
ing are recommended, regarding the execution of the modified 
CHAMP ION/ 2/E/FIX program. 

a. The number of TDMA traverses should be reduced from 
five to three. 

b. The number of iterations may remain the same or be 
reduced to one half without significant reduction in the 
accuracy of predictions. The use of any number of iterations 
greater than 200 is not expected to promote the accuracy 
significantly. 

c. While establishing the grid, the nodes which are 

closer to the solid fuel grain wall should be selected in 

* 

such a way as to result in values greater than 11.5 (at 
least for most of them) . This will increase the accuracy 
of predictions and will enable the investigator to use fewer 
grid lines in both the x and y directions. 
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APPENDIX A 



C/ASPI#5 JOB (3238,1085) , 'GENMIX' ,CLASS=C CHAOO 

C/*MAIN LINES=(12) CHAOO 

C/ EXEC FORTXCG,REGION=1024K CHAOO 

C/FORT. SYSPRINl DD DUMMY CHAOO 

C/FORT. SYS IN DD * CHAOO 

C ############################ CHAOO 

C ANDREA MH SVYSHS AYTO TO FAKELLO TON XREIAZOMAI GIA 

TH THESIS CHAOO 

C ################################ CHAOO 

BLOCK DATA CHAOO 

COMMON CHAOO 

1/CASE51/RIN, lYIN, lYINPl , IXIN, IXINMl ,TAULW (42) , 

XPUSLW (42), CTAULW , CHAO 0 

1 CXPLW CHAOO 

1/CAS E52/1XIN2 , IYIN2 , NXRl , NYRl , NXRlMl , NYRlMl , 

NYRIM2, CHAOO 

1 CYPTW1,CTAUT1,CXPLW2,CTAUL2,HWALL2 CHAOO 

1/CASE 1/UINLET , FLOWIN , RP IPE , XPIPE , FXSTEP , HINLET , 

HWALL CHAOO 

2/DNY/ DYG(42) ,CYV(42) ,FV(42) ,FVNODE(42) ,R(42) , 

RDYG(42) ,RDYV(42) CHAOO 

2, RSYG(42) ,RSYV(42) ,RV(42) ,RVCB(.42) ,RVSQ(42) , 

SYG(42) ,SYGTR(42) . CHAOO 

2,SYV(42) ,SYVTRV(42) ,Y (42)‘,YV(42) CHAOO 

3/DNYONX/AE(42) ,AN(42) ,AP (42) ,AS (42) ,AW(42) ,C (42) , 

D(42) ,DIFE(42) CHAOO 

3, DIFN(42) ,0017(42) ,DIFW(42) ,DU(42) ,DV(42) ,EMUE(42) , 

EMUN(42) CHAOO 

3,EMUW(42) ,USTAR(42) ,VSTAR(42) CHAOO 

3,PHIOLD(,42) ,RHCE(42) ,RHON(42) ,RHOW(42) ,S (42) , 

SPRIME(42) CHAOO 

3, VOLUME (42) , CONN (42) , CONE (42) ,CONW(42) CHAOO 

4/DNX/DXG(47) ,DXU(47) ,FAREA(47) ,FU(47) CHAOO 

4, FUNODE(.47) ,KOUNT(47) ,RDXG(47i CHAOO 

4, RDXU(47) ,RSXG(47) ,RSXU(47) , STORE (.47) ,SXG(.47) , 

SXU(47) ,X(47) ,XU(47) CHAOO 

COMMON CHAOO 

5/DJPHI/ lEW(IC) ,ILAST(10) , IKON (.10) ,IXNY(.10) , 

I ZERO (.10) CHAOO 

5, KS0LVE(11) ,KRS(10) , RELAX (10) ,RSREF(10) ,RSSUM(10) , 

TITLE (10) CHAOO 

COMMON CHAOO 

6/DO/CCHECK, DP , FLOWPC , FLOWST , FLOWUP , GREAT , IPLRS , 

IPREF,IPRINT CHAOO 

6 , IWEEEP , IX , IXINY , IXlNYl , IX2NY2 , IXMON , IXPl , IXPREF , 

IXU,IYMON,IYPREF CHAOO 
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6 , JEMU , JH , JLAST , JP , JPP , JRHO , JU , JV , JVPl , KALLED , 

KINPRI^KEHOMU CKAOO 

6 , KTEST , LABPHI , LSWEEP , NSOLVE , NTRAV , KRAD CHAO 0 

6 , lYF , lYFl , lYL , lYL^ll , lYLPl ,NTDMA , NUMCOL CHAO 0 

6,N) ,NXMA) ,NXM1,NXM2 ,NXYG,NXYP,NXYU,NXYV CHAOO 

6,NYNYKA) ,NYM1,NYM2 ,PI,RSCHEK,RSMAX,TINY CHAOO 

6,IPFRS1,IPFRS2,KIN,KEX CHAOO 

COMMON/PROP , EMUREF , PREFF (10), PRL (10), PRT (10), 

RHOREF CHAOO 

COMMON/D2D1/ARSL (47,10) ,RSLINE (47,10) CHAOO 

COMMON/D2D2/U(1560) ,V(1560) ,TKEC1600) , TED (1600) , 

H(1600) ,PP(40) , CHAOO 

7P(1444) ,RHO(1600) ,EML(1600) ,RSF(40,40) ,T(1600) CHAOO 

DIMENSION F (14200) CHAOO 

COMMON CHAOO 

9/TURB/Cl , 02 , CD , SCRTCD , CD2 5 , ECONST , CTAUTW , CYPTW , 

TAUTW (47), GENK (.42), CHAO 0 

9 YPUSTW (47), FACTKE , FACTED , JTKE , JTED , CAP PA , STANTl'7 (47), 

HBULKC47), CHAOO 

9 RUN , CP , GC , HINR , HWALLR , PREFR , BARM , JT , HOHECK CHAO 0 

9/COMB/AWALL (42) ,FLUX (42), FLOWAL (4 2 ) , BP (.4 2 ) , 

RR (42) ,TEMP1 (.42) , CHAOO 

9EMUF , RHCF , HUF , HTR , TREF , STC , HC ,FM2 IN ,FOXIN , CHAOO 

9FMCLFU , FMCLPR , FMCLOX , FMOLN2 CHAO 0 

9/CHEMCO/FFU(42) ,FOX(42) ,FN2 (.42) , CHAOO 

9FLFUIT , FLFLOT , FLOXIT , FLOXCT , FLN2 IT , FLN2 ) T , CHAO 0 

9PCEFU,PCECX,PCEN2,RRAVE,LINE CHAOO 

REAL*8 TL(12) CHAOO 

LOGICAL*! LTG(3) CHAOO 

COMMON ‘ CHAOO 

D/OPLOT/TL CHAOO 

D, KPLOT(ll) ,NINT(11) ,PHI (40,40) ,CL(25) ,XP(47) , 

SCALEW,SCALEH,LTG CHAOO 

COMMON/RAD/TWALL , TFLAME , S I GMA , EPS IW , IRAD CHAOO 

E, IRADS, RADI, FUEL (4 0) ,QRAD(40) ,QCON(4 0) ,F^!ASF (.40). , 

VOLL(40) CHAOO 

DIMENSION DIFS (42) ,EMUS (4 2) ,RHOS (.42) ,CONS (42) CHAOO 

EQUIVALENCE (CONS (2) ,CONN (1) ) CHAOO 

EQUIVALENCE (DIFS (2) ,DIFN (1) ) , (EMUS (2) ,EMUN (1) ) CHAOO 
C EQUIVALENCE (RHOS(2) ,RHON(-l) ) , (AREAE, AREAW) CHAOO 

EQUI VALNECE ( RHOS ( 2 ) , RHON (.1 ) ) CHAO 0 

EQUIVALENCE (F(1),U(1)) CHAOO 

DIMENSION A (42) ,B (42) CHAOO 

EQUIVALENCE (A CD ,AN ( 1) ) , (B (1) , AS (1) ) CHAOO 

CHAPTER 1 GENERAL, FLOW PARAMETERS CHAOO 

DATA GREAT, TINY, PI/1. E30, l.E-30, 3.1415926/ CHAOO 

C CHAOO 

DATA RPIPE,XPIPE,UINLET,HINLET,HWALL/ CHAOO 

10.026924,0.4628,92.220,0.0,1.0/ CHAOO 

CHAOO 
CHAOO 

DATA HWALL 2/1.0/ CHAOO 

DATA IXIN,IYIN/2,7/ CHAOO 
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DATA IXIN,IYIN/2, 7/ CHAOO 

C DATA KTEST/0/ CHAOO 

DATA KTEST/1/ CHAOO 

CHAPTER 2 GRID CHAOO 

DATA NXMAX,NYMAX/40,40/ CHAOO 

C CHAOO 

DATA NX, NY/ 4 0,3 3/ CHAOO 

DATA KRAD/2/ CHAOO 

C KRAD=2 MEANS CYLINDRICAL POLAR COORDINATES. CHAOO 

C IF ONE WISHES THE PROGRAM TO CALCLXATE THE GRID, CHAOO 

C SET X(.l) OR Y(l) NOT EQUAL TO 0 . 0 CHAOO 

C CHAOO 

C CHAOO 



DATA X/0.0, .01, .023, .038, .055, .07333, .09166, 

.10999, .12832, CHAOO 

A. 14665, .16498, .18331, .20164, .21997, . 2383 ,. 25663 , 

.27343, .28468 CHAOO 

B. 29218, .29718, .30055, .3028, .3043, .3053, CHAOO 

C. 3066, .3081, .3096, .31185, .31522, . 32022 ,. 32772 , 

.33892, .35392, CHAOO 

D. 36892, .38292, .39892, . 4139 2 , . 4289 2 , .44192, 

.4543, .4643, CHAOO 

EO.,0. ,0.,0. ,0.,0./ CHAOO 

DATA Y/0.0, .00061, .00141, .00271, .00401, .00531, 

.00631, .00731, CHAOO 

A. 00861, .01011, .012235, .01436, . 01586 ,. 01686 , 

.01753, .01798, .01828, CHAOO 

B. 01853, .01873, .01888, . 019 00 ,. 01910 , .01925, CHAOO 

C. 01955, .02015, .02112, .02215, . 02315 ,. 02415 ,. 02515 , CHAOl 

D. 02565, .026260, .026924,0. ,0.,0.,0.,0.,0.,0.,0.,0./ CHAOl 

C CHAOl 

DATA FXSTEp/1.5/ CHAOl 

CHAPTER 3 VARIABLES CHAOl 

C VARIABLES SOLVED BY TDMA CHAOl 

DATA JU,JV,JTKE,JTED,JH,JPP/1,2,3,4,5,6/ CHAOl 

C AUXILIARY VARIABLES CHAOl 

DATA JP, JRHO,JEMU,JT,JLAST/7,8,9 ,10,10/ CHAOl 

C the number of DEPENDENT VARIABLES SOLVED CHAOl 

DATA NSOLVE/5/ CHAOl 

DATA KSOLVE/9*!, 2*0/ CHAOl 

DATA KRS/6*1,4*0/ CHAOl 

C CHAOl 

DATA RELAX/0.3,0.3,0.5,0.5,1.,1.,0.1.,01,0.5,1./ CHAOl 
DATA TITLE/10 *4 HXXXX/ CHAOl 

CHAPTER 4 PROPERTY DATA CHAOl 

DATA RHOREF,EMUREF/1.403,3.034E-5/ CHAOl 

DATA PRL,PRT/12*1. , .9, 7*1./ CHAOl 

DATA PREFF/3*1.0, 1.3, 6*1.0/ CHAOl 

DATA C1,C2,CD,CAPPA,ECONST/1.43,1.92,0.09,0.4 ,9 .0/ CHAOl 
C EPSIW IS WALL EMISSIVITY,RAD1 IS EMP . CONST .FOR CHAOl 

C EMISSGAS/UNBURNED FUEL ON A LINE CHAOl 

DATA SIGMA, EPS IW,TWALL,TFLAME, RADI/ CHAOl 
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15.669E-08,0.85,600. , 319 0 . , 3 . 0000/ 

C IRAD=0,BP NOT MODIFIED FOR RADIATION 

C IR/D=1 WITH RADIATION 

DATA IRAD/0/ 

C DATA IRAD/1/ 

CHAPTER 5 STARTING PREPARATIONS 

DATA IXPREF , lYPREF/ 2 , 2 / 

DATA IXMON,IYMON/22,24/ 

DATA KINPRI/0/ 

DATA F,DV/14200*0 .0 ,42*0.0/ 

DATA FACTKE,FACTED/0 .005,0 .03/ 

CHAPTER 6 STEP CONTROL 

CHAPTER 7 BOUNDARY CONDITIONS 

DATA KIN,KEX/3,1/ 

CHAPTER 8 ADVANCE 

DATA NTDMA/5/ 

DATA RSCHEK/0.001/ 

CHAPTER 9 COMPLETE 

C 

C SUBROUTINE DENS IT CALCULATES VARIABLE DENSITY 

ON EACH LINE. 

C RUN=UNIVERSAL GAS CONSTANT 

C CP=C SUB P 

C PREFR IS THE REFERENCE PRES ADDED TO P(I) TO 

CALCULATE P ACTUAL 
C 

DATA FLOWAL/42*0 .0/ 

DATA RUN, OP, GC/8314 .0,1004 .0 ,1.0/ 

C 

' DATA HINP, HWALUR, PREFR, BARM/2. 6570E06, 4 .826E0 5, 
3.152E05,28.964/ 

DATA HTR,FOXIN,FN2IN,STC,HC/2.300E16,0 .233, 
0.767, 1.92, 2. 443E07/ 

DATA FMOLFU , FMCLPR , FMOLOX , FMOLN2/1 0 0.0,28.97, 
32.0,28.0/ 

DATA RHOF,EMUF/1199 . ,2 .43E-05/ 

DATA TREF/600.0/ 

C 

C SET VARIABLE 'LINE' = IX ON THE LINE ON WHICH 

C CONSERVATION CALCULATIONS DESIRED. 

DATA LINE/29/ 

C 

CHAPTER 10 ADJUST 

CHAPTER 11 PRINT 

DATA NUMCOL/10/ 

DATA IPLRS,IPRINT/20 0,200/ 

C ### CORRECTION ON BASIC PROGRAM 
DATA IPLRS,IPRINT/200,200/ 

CHAPTER 12 DECIDE 

DATA CCHECK/0 .001/ 

DATA LSWEEP.200/ 



CHAOl 

CHAOl 

CHAOl 

CHAOl 

CHAOl 
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CHAOl 

CHAOl 

CHAOl 
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CHAOl 
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o o o o o o on o o 



C######################################## CHAOl 

C DATA LSWEEP/2/ CHAOl 

C############################################## CHAOl 

C CHAOl 

END CHAOl 

C PROGRAM MAIN ( OUTPUT, TAPE 6 =OUTPUT) CHAOl 

C CCHaOl 

C FOR CHAMPION 2/E/FIX, CASE50, (HTS762 8, C50) CHAOl 

C IDENT CASE 51 CHAOl 

C 29.7,77 CHAOl 

C CHAOl 

C CCHAOl 

C CCHAO 1 

C/JULY 1977/ CCHAOl 

CCHAO 1 

A DETAILED DESCRIPTION OF THE PROGRAM CAN BE CCHAO 1 

FOUND IN THE REPORT A GENERAL COMPUTER CCHAOl 

PROGRAM FOR TWO-DIMENSIONAL ELLIPTIC FLOWS CCHAOl 

BY W.M. PUN AND D.B. SPALDING, HTS/76/2, CCHAOl 

MECHANICAL ENGINEERING DEPARTMENT, IMPERIAL CCHAO 1 
COLLEGE, LONDON CCHAOl 



CCHAO 1 
CCHAO 1 
CCHAOl 

CHAPTER 000000000 DECLARATIONS 0 0 0 v CHAOl 



C CHAOl 

COMMON CHAOl 

1/CASE51/RIN , lYIN , lYINPl , IXIN , IXINml , TAULW (42), 

XPUSLW(42) ,CTAULW, CHAOl 

1 CXPLW CHAOl 

1/CASE52/IXIN2 , IYIN2 , NXRl , NYRl , NXRlMl , NYRlMl , 

NYR1M2 , CHAO 1 

1 CYPTVa , CTAUTl , CXPLW2 , CTAUL2 , HWALL2 CHAO 1 

1/CASEl/UINLET , FLOWIN , RPIPE ,XPIPE , FXSTEP , HINLET , 

HWALL CHAOl 

2/DNY/ DYG(42) ,DYV(42) ,FVC42) ,FVNODE(42) ,R(42I , 

RDYGC42) ,RDYV(42) CHA02 

2, RSYG(42) ,RSYV(42) ,RV(42) ,RVCB(42) ,RVSQ(42) , 

SYGC42) ,SYGTR(42) CHA02 

2,SYV(42) ,SYVTRV(42) ,Y(42) ,YV(.42) CHA02 

3/DNYONX/AE(42) ,AN(42) ,AP(42) ,AS(42), ,AW(42) ,C(42)_ , 

D(42) ,DIFE(42) CHA02 

3, DIFN(42) ,DUW(42) ,DIFW(42)_ ,DUC42), ,DV(.421 ,EMUEC421 , 

EMUNC42) CHAO 2 

3,EMUWC42) ,USTAR(42) ,VSTAR(,42) CHA02 

3,PHIOLD(.42) ,RH0EC42) ,RH0NC4 2I ,RH0WC42) ,SC42) , 

SPRIMEC42) CHAO 2 

3, VOLUME(42) , CONN (42) , CONE (42) ,CONW(42) CHA02 

4/DNX/DXG (.47) , DXU (.4 7) , FAREA (47) , FU (.4 7 ) CHAO 2 

4, FUNODE(47) ,KOUNT(47) ,RDXG(47) CHA02 

4,RDXU(47) ,RSXG(47) ,RSXU(47) ,STORE(47) ,SXG(47). , 

SXU(47) ,X(47) ,XU(47) CHA02 
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CHAO 2 



5/DJPHI/ lEW(lO) ,ILAST(10) ,IMON(10) ,IXNY(10) , 

I ZERO (10) 

5, KS0LVE(11) ,KRS (10) , RELAX (10) ,RSREF(10) ,RSSUM(10) , 

TITLE (10) CHAO 2 

COMMON CHAO 2 

6/DO/CCHECK , DP , FLOWPC , FLOWST , FLOWUP , GREAT , IPLRS , 

IPREF,IPRINT CHAO 2 

6 , ISWEEP , IX , IXINY , IXlNYl , IX2NY2 , IXMON , IXPl , IXPREF , 

IXU , I YMON , I YPREF CHAO 2 

6 , JEMU , JH , JLAST , JP , JPP , JRH 0 , JU , JV , JVPl , KALLED , 

KINPRI , KRHOMU CHAO 2 
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